
rm( list = ls() )
library(MatchIt)
n_units <- 1000
n_covars <- 5
n_periods <- 2 


coef_vec <- matching2_coef_vec <- matching_coef_vec <- c()

for(i in 1:1000){ 
print(i) 
    
  X <- matrix(rnorm(n_units * 2 * (n_covars + 1)), 
              nrow = n_units * 2, 
              ncol = n_covars)
  treatment_vec <- rep(rbinom(n_units, 
                          size = 1, 
                          prob = 0.5), n = n_periods)
  post_vec <- rep(0, times = n_units * n_periods)
  post_vec[-c(1:n_units)] <- 1
  
  interaction_vec <- post_vec * treatment_vec
  treat_coef <- 1
  post_coef <- 1
  interaction_coef <- 2
  x_coefs <- as.matrix(  rep(0.25, times = n_covars) ) 
  sigma2 <- 0.25
  Y_obs <- treat_coef * treatment_vec + 
              post_coef * post_vec + 
                interaction_coef * interaction_vec+ 
                      X %*% x_coefs + 
                        rnorm(2 * n_units, mean = 0, sd = sqrt(  sigma2)  ) 
  Y_obs <- c( Y_obs  ) 
  full_df <- as.data.frame( cbind(Y_obs, treatment_vec, interaction_vec, post_vec, X) ) 
  full_df$ID <- rep(1:n_units, times = 2)
  
  #without matching on Y_obs (pre)
  coef_vec[i] <- coef( summary(  lm(Y_obs~interaction_vec + 
                                               treatment_vec + post_vec + 
                                               V5 + V6 + V7 + V8 + V9, 
                                             data = full_df) )  )[2,1]
  
  
  #with matching on Y_obs (pre)
  wide_df <- tapply(1:nrow(full_df), full_df$ID, function(x){
          Contrib1 <- full_df[x[1],]
          names(Contrib1) <- paste("Time1", colnames(Contrib1), sep = "_")
          
          Contrib2 <- full_df[x[2],]
          names(Contrib2) <- paste("Time2", colnames(Contrib2), sep = "_")
          return(cbind(Contrib1, Contrib2))
          })
  wide_df <- do.call(rbind, wide_df)
  match_list <- matchit(Time2_treatment_vec~Time1_Y_obs, 
                        data = wide_df, 
                        distance = "logit", 
                        ratio = 3, 
                        replace = T, 
                        caliper = 0.2)
  match_wts <- rep(match_list$weights, 2)
  matching_coef_vec[i] <- coef( summary(  lm(Y_obs~interaction_vec + 
                                               treatment_vec + post_vec + 
                                               V5 + V6 + V7 + V8 + V9,
                                             weights = match_wts, 
                                             data = full_df) )  )[2,1]
  if(T == T){ 
    keep_IDs <- wide_df$Time2_ID [match_list$weights > 0]
    full_df_reduced <- full_df[full_df$ID %in% keep_IDs, ]
    matching_coef_vec[i] <- coef( summary(  lm(Y_obs~interaction_vec + 
                                               treatment_vec + post_vec + 
                                               V5 + V6 + V7 + V8 + V9,
                                             data = full_df_reduced) )  )[2,1]
  }
  
  #matching 2 
  match_list2 <- matchit(Time2_treatment_vec~Time1_V5 + Time1_V6 + Time1_V7 + Time1_V8 + Time1_V9, 
                        data = wide_df, 
                        distance = "logit", 
                        ratio = 3, 
                        replace = T, 
                        caliper = 0.2)
  match_wts2 <- rep(match_list2$weights, 2)
  matching2_coef_vec[i] <- coef( summary(  lm(Y_obs~interaction_vec + 
                                               treatment_vec + post_vec + 
                                               V5 + V6 + V7 + V8 + V9,
                                             weights = match_wts2, 
                                             data = full_df) )  )[2,1]
  
  if(T == T){ 
    keep_IDs2 <- wide_df$Time2_ID [match_list2$weights > 0]
    full_df_reduced2 <- full_df[full_df$ID %in% keep_IDs2, ]
    matching2_coef_vec[i] <- coef( summary(  lm(Y_obs~interaction_vec + 
                                                 treatment_vec + post_vec + 
                                                 V5 + V6 + V7 + V8 + V9,
                                               data = full_df_reduced2) )  )[2,1]
  }
} 



pdf("~/Dropbox/Collaborations/Traffic/MonteCarlo_TechnicalAppendix.pdf", 
    width = 7, height = 8)
par(mfrow = c(2,1))
cex_main <- 2 
cex_lab <- 1.5
my_xlim <- c(min(coef_vec, matching_coef_vec, matching2_coef_vec), 
             max(coef_vec, matching_coef_vec, matching2_coef_vec) )
plot(density(matching2_coef_vec), 
     xlim = my_xlim, 
     cex.main = cex_main, 
     cex.lab = cex_lab, 
     xlab = "Estimate Across 1000 Monte Carlo Simulations", 
     main = "Matching on Pre-Treatment Covariates")
abline(v = interaction_coef, 
       lty = 90, 
       lwd = 3, 
       col = "red")
text(x = 1.965, 
     y = 2, 
     labels = "Truth")
arrows(x0 = 1.98, y0 = 2, x1 = 2, y1 = 2.5, length = 0.15, lty = 1)

plot( density(matching_coef_vec), 
      xlim = my_xlim, 
      cex.main = cex_main, 
      cex.lab = cex_lab, 
      xlab = "Estimate Across 1000 Monte Carlo Simulations", 
      main = "Matching on Pre-Intervention Outcome")
abline(v = interaction_coef, 
       lty = 90, 
       lwd = 3, 
       col = "red")
text(x = 1.965, 
     y = 2, 
     labels = "Truth")
arrows(x0 = 1.98, y0 = 2, x1 = 2, y1 = 2.5, length = 0.15, lty = 1)
dev.off()

